LAMMPS (10 Sep 2025 - Development - patch_10Sep2025-520-gd2199aa57b-modified)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread.
  using 1 OpenMP thread(s) per MPI task
# Analytical calculation
# of Born matrix

# Note that because of cubic symmetry and central forces, we have:
# C11, pure axial == positive mean value: 1,2,3
# C44==C23, pure shear == positive mean value, exactly match in pairs: (4,12),(5,8),(6,7)
# C14==C56, shear/axial(normal) == zero mean, exactly match in pairs: (9,21),(14,20),(18,19)
# C15, shear/axial(in-plane) == zero mean: 10,11,13,15,16,17

# adjustable parameters

units           real
variable        nsteps index 10000     # length of run
variable        nthermo index  1000    # thermo output interval
variable        nlat equal 5           # size of box
variable        T    equal 60.         # Temperature in K
variable        rho  equal 5.405       # Lattice spacing in A

atom_style      atomic

lattice         fcc ${rho}
lattice         fcc 5.405
Lattice spacing in x,y,z = 5.405 5.405 5.405
region          box block 0 ${nlat} 0 ${nlat} 0 ${nlat}
region          box block 0 5 0 ${nlat} 0 ${nlat}
region          box block 0 5 0 5 0 ${nlat}
region          box block 0 5 0 5 0 5
create_box      1 box
Created orthogonal box = (0 0 0) to (27.025 27.025 27.025)
  1 by 2 by 2 MPI processor grid
create_atoms    1 box
Created 500 atoms
  using lattice units in orthogonal box = (0 0 0) to (27.025 27.025 27.025)
  create_atoms CPU = 0.000 seconds

mass            * 39.948

velocity        all create ${T} 87287 loop geom
velocity        all create 60 87287 loop geom
velocity        all zero linear

pair_style      lj/cut 12.0
pair_coeff      1 1 0.238067 3.405

neighbor        0.0 bin
neigh_modify    every 1 delay 0 check no

variable vol equal vol
thermo 100
fix aL all ave/time 1 1 1 v_vol ave running
fix NPT all npt temp $T $T 100 aniso 1. 1. 1000 fixedpoint 0. 0. 0.
fix NPT all npt temp 60 $T 100 aniso 1. 1. 1000 fixedpoint 0. 0. 0.
fix NPT all npt temp 60 60 100 aniso 1. 1. 1000 fixedpoint 0. 0. 0.

run 2000
Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
  update: every = 1 steps, delay = 0 steps, check = no
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 12
  ghost atom cutoff = 12
  binsize = 6, bins = 5 5 5
  1 neighbor lists, perpetual/occasional/extra = 1 0 0
  (1) pair lj/cut, perpetual
      attributes: half, newton on
      pair build: half/bin/atomonly/newton
      stencil: half/bin/3d
      bin: standard
Per MPI rank memory allocation (min/avg/max) = 2.649 | 2.649 | 2.649 Mbytes
   Step          Temp          E_pair         E_mol          TotEng         Press          Volume    
         0   60            -974.97988      0             -885.73443     -1549.5634      19737.726    
       100   33.504753     -943.70789      0             -893.87211      968.3116       18369.007    
       200   26.32504      -930.72444      0             -891.56794      314.3869       18935.888    
       300   41.132299     -931.36826      0             -870.18708     -617.04719      19632.688    
       400   41.183387     -912.49497      0             -851.2378       485.3659       19105.509    
       500   57.063649     -914.95616      0             -830.07831      10.213449      19451.639    
       600   57.720228     -897.82416      0             -811.96969     -122.66415      19794.9      
       700   55.117494     -878.27174      0             -796.28864      169.40743      19853.365    
       800   65.753909     -877.79094      0             -779.98698     -140.35502      20145.428    
       900   61.5708       -870.58334      0             -779.00144     -6.2782272      20132.09     
      1000   59.71058      -873.44425      0             -784.62929      18.589753      20063.414    
      1100   59.887921     -878.67927      0             -789.60053     -53.049476      20034.152    
      1200   59.121142     -881.97593      0             -794.03771      31.110495      19909.638    
      1300   60.931065     -887.19784      0             -796.5675       6.5972135      19859.309    
      1400   60.086091     -890.68151      0             -801.308       -26.605388      19827.013    
      1500   58.240958     -892.81245      0             -806.18344      32.781814      19743.494    
      1600   59.997061     -896.10853      0             -806.86745      1.8694388      19731.408    
      1700   62.500416     -900.56177      0             -807.59714     -58.405946      19717.51     
      1800   60.950957     -905.28688      0             -814.62695      13.287333      19580.499    
      1900   57.643425     -908.35702      0             -822.61679      80.300088      19476.931    
      2000   54.860033     -904.52282      0             -822.92268      28.955012      19564.956    
Loop time of 1.30444 on 4 procs for 2000 steps with 500 atoms

Performance: 132.471 ns/day, 0.181 hours/ns, 1533.228 timesteps/s, 766.614 katom-step/s
99.8% CPU use with 4 MPI tasks x 1 OpenMP threads

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 0.43945    | 0.4508     | 0.46521    |   1.5 | 34.56
Neigh   | 0.65198    | 0.66196    | 0.6713     |   1.1 | 50.75
Comm    | 0.13772    | 0.16144    | 0.18284    |   4.0 | 12.38
Output  | 0.00025968 | 0.00032781 | 0.00050556 |   0.0 |  0.03
Modify  | 0.023953   | 0.024197   | 0.024365   |   0.1 |  1.85
Other   |            | 0.005718   |            |       |  0.44

Nlocal:            125 ave         129 max         121 min
Histogram: 1 0 0 1 0 0 1 0 0 1
Nghost:           1737 ave        1741 max        1733 min
Histogram: 1 0 0 1 0 0 1 0 0 1
Neighs:        11600.2 ave       12136 max       11091 min
Histogram: 1 0 0 1 0 1 0 0 0 1

Total # of neighbors = 46401
Ave neighs/atom = 92.802
Neighbor list builds = 2000
Dangerous builds not checked

unfix NPT

variable newL equal "f_aL^(1./3.)"
change_box all x final 0 ${newL} y final 0. ${newL} z final 0. ${newL} remap units box
change_box all x final 0 26.9850129540338 y final 0. ${newL} z final 0. ${newL} remap units box
change_box all x final 0 26.9850129540338 y final 0. 26.9850129540338 z final 0. ${newL} remap units box
change_box all x final 0 26.9850129540338 y final 0. 26.9850129540338 z final 0. 26.9850129540338 remap units box
Changing box ...
  orthogonal box = (0 0 0) to (26.985013 26.94049 27.008247)
  orthogonal box = (0 0 0) to (26.985013 26.985013 27.008247)
  orthogonal box = (0 0 0) to (26.985013 26.985013 26.985013)

unfix aL

reset_timestep 0

# Conversion variables
variable kb        equal 1.38065e-23    # J/K
variable Myvol     equal "vol*10^-30" # Volume in m^3
variable kbt       equal "v_kb*v_T"
variable Nat       equal atoms
variable Rhokbt    equal "v_kbt*v_Nat/v_Myvol"
variable at2Pa     equal 101325
variable kcalmol2J equal "4183.9954/(6.022e23)"
variable C1        equal "v_kcalmol2J/v_Myvol" # Convert Cb from energy to pressure units
variable C2        equal "v_Myvol/v_kbt"       # Factor for Cfl terms
variable Pa2GPa    equal 1e-9

# Born compute giving <C^b> terms
compute     born all born/matrix
# The six virial stress component to compute <C^fl>
compute     VIR all pressure NULL virial
variable s1 equal "-c_VIR[1]*v_at2Pa"
variable s2 equal "-c_VIR[2]*v_at2Pa"
variable s3 equal "-c_VIR[3]*v_at2Pa"
variable s6 equal "-c_VIR[4]*v_at2Pa"
variable s5 equal "-c_VIR[5]*v_at2Pa"
variable s4 equal "-c_VIR[6]*v_at2Pa"
variable press equal press


# Average of Born term and vector to store stress
# for post processing
fix CB all ave/time 1 ${nthermo} ${nthermo} c_born[*] ave running file born.out overwrite
fix CB all ave/time 1 1000 ${nthermo} c_born[*] ave running file born.out overwrite
fix CB all ave/time 1 1000 1000 c_born[*] ave running file born.out overwrite
fix CPR all ave/time 1 1 1 c_VIR[*] file vir.out
fix APR all ave/time 1 1 1 v_press ave running
fix VEC all vector 1 v_s1 v_s2 v_s3 v_s4 v_s5 v_s6

thermo      ${nthermo}
thermo      1000
thermo_style    custom step temp press f_APR c_born[1] f_CB[1] c_born[12] f_CB[12] c_born[4] f_CB[4]
thermo_modify line multi

fix     1 all nvt temp $T $T 100
fix     1 all nvt temp 60 $T 100
fix     1 all nvt temp 60 60 100

run         ${nsteps}
run         10000
Generated 0 of 0 mixed pair_coeff terms from geometric mixing rule
Neighbor list info ...
  update: every = 1 steps, delay = 0 steps, check = no
  max neighbors/atom: 2000, page size: 100000
  master list distance cutoff = 12
  ghost atom cutoff = 12
  binsize = 6, bins = 5 5 5
  2 neighbor lists, perpetual/occasional/extra = 1 1 0
  (1) pair lj/cut, perpetual
      attributes: half, newton on
      pair build: half/bin/atomonly/newton
      stencil: half/bin/3d
      bin: standard
  (2) compute born/matrix, occasional
      attributes: full, newton on
      pair build: full/bin/atomonly
      stencil: full/bin/3d
      bin: standard
Per MPI rank memory allocation (min/avg/max) = 3.28 | 3.28 | 3.28 Mbytes
------------ Step              0 ----- CPU =            0 (sec) -------------
Step     =              0 Temp     =        54.8600 Press    =       -74.9874 
f_APR    =       -74.9874 c_born[1] =      9269.2848 f_CB[1]  =         0.0000 
c_born[12] =      5315.9657 f_CB[12] =         0.0000 c_born[4] =      5315.9657 
f_CB[4]  =         0.0000
------------ Step           1000 ----- CPU =     2.682027 (sec) -------------
Step     =           1000 Temp     =        56.1005 Press    =        93.2723 
f_APR    =       104.8993 c_born[1] =      9899.7837 f_CB[1]  =      9738.6101 
c_born[12] =      5515.1091 f_CB[12] =      5497.0605 c_born[4] =      5515.1091 
f_CB[4]  =      5497.0605
------------ Step           2000 ----- CPU =     5.372834 (sec) -------------
Step     =           2000 Temp     =        59.2245 Press    =       164.2019 
f_APR    =        90.3569 c_born[1] =     10127.7921 f_CB[1]  =      9762.5811 
c_born[12] =      5582.5570 f_CB[12] =      5512.1775 c_born[4] =      5582.5570 
f_CB[4]  =      5512.1775
------------ Step           3000 ----- CPU =     8.024821 (sec) -------------
Step     =           3000 Temp     =        62.3421 Press    =       138.9920 
f_APR    =       103.8105 c_born[1] =      9784.0632 f_CB[1]  =      9805.4617 
c_born[12] =      5558.8869 f_CB[12] =      5514.3234 c_born[4] =      5558.8869 
f_CB[4]  =      5514.3234
------------ Step           4000 ----- CPU =     10.64483 (sec) -------------
Step     =           4000 Temp     =        58.1487 Press    =       136.4624 
f_APR    =       107.3823 c_born[1] =      9996.9126 f_CB[1]  =      9821.2552 
c_born[12] =      5439.8328 f_CB[12] =      5518.8966 c_born[4] =      5439.8328 
f_CB[4]  =      5518.8966
------------ Step           5000 ----- CPU =     13.34843 (sec) -------------
Step     =           5000 Temp     =        60.3351 Press    =       -21.9127 
f_APR    =       101.0136 c_born[1] =      9631.2101 f_CB[1]  =      9798.9441 
c_born[12] =      5352.0607 f_CB[12] =      5519.0345 c_born[4] =      5352.0607 
f_CB[4]  =      5519.0345
------------ Step           6000 ----- CPU =     16.11892 (sec) -------------
Step     =           6000 Temp     =        57.7282 Press    =       157.6900 
f_APR    =       100.3234 c_born[1] =     10030.8342 f_CB[1]  =      9800.1967 
c_born[12] =      5740.8579 f_CB[12] =      5521.0716 c_born[4] =      5740.8579 
f_CB[4]  =      5521.0716
------------ Step           7000 ----- CPU =     18.80235 (sec) -------------
Step     =           7000 Temp     =        59.5948 Press    =       135.9651 
f_APR    =       101.9810 c_born[1] =     10058.4542 f_CB[1]  =      9817.6242 
c_born[12] =      5511.2599 f_CB[12] =      5521.8274 c_born[4] =      5511.2599 
f_CB[4]  =      5521.8274
------------ Step           8000 ----- CPU =     21.45934 (sec) -------------
Step     =           8000 Temp     =        61.7547 Press    =        41.1860 
f_APR    =       105.3152 c_born[1] =      9554.9451 f_CB[1]  =      9830.6704 
c_born[12] =      5292.0749 f_CB[12] =      5520.3055 c_born[4] =      5292.0749 
f_CB[4]  =      5520.3055
------------ Step           9000 ----- CPU =     24.13119 (sec) -------------
Step     =           9000 Temp     =        59.2478 Press    =       146.6808 
f_APR    =       104.7030 c_born[1] =      9860.2050 f_CB[1]  =      9823.1198 
c_born[12] =      5543.4620 f_CB[12] =      5518.4889 c_born[4] =      5543.4620 
f_CB[4]  =      5518.4889
------------ Step          10000 ----- CPU =     26.77791 (sec) -------------
Step     =          10000 Temp     =        59.4645 Press    =        96.1963 
f_APR    =       107.2322 c_born[1] =      9779.2285 f_CB[1]  =      9836.0535 
c_born[12] =      5617.0229 f_CB[12] =      5525.6584 c_born[4] =      5617.0229 
f_CB[4]  =      5525.6584
Loop time of 26.7779 on 4 procs for 10000 steps with 500 atoms

Performance: 32.265 ns/day, 0.744 hours/ns, 373.442 timesteps/s, 186.721 katom-step/s
99.7% CPU use with 4 MPI tasks x 1 OpenMP threads

MPI task timing breakdown:
Section |  min time  |  avg time  |  max time  |%varavg| %total
---------------------------------------------------------------
Pair    | 2.0638     | 2.1385     | 2.1975     |   3.3 |  7.99
Neigh   | 3.1134     | 3.1287     | 3.1565     |   1.0 | 11.68
Comm    | 0.67685    | 0.76518    | 0.8541     |   7.2 |  2.86
Output  | 0.012628   | 0.012663   | 0.012763   |   0.1 |  0.05
Modify  | 20.711     | 20.712     | 20.714     |   0.0 | 77.35
Other   |            | 0.02053    |            |       |  0.08

Nlocal:            125 ave         131 max         119 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Nghost:           1737 ave        1743 max        1731 min
Histogram: 2 0 0 0 0 0 0 0 0 2
Neighs:        11548.8 ave       12134 max       10971 min
Histogram: 1 1 0 0 0 0 0 0 1 1
FullNghs:      23097.5 ave       24209 max       21987 min
Histogram: 2 0 0 0 0 0 0 0 0 2

Total # of neighbors = 92390
Ave neighs/atom = 184.78
Neighbor list builds = 10000
Dangerous builds not checked

# Compute vector averages
# Note the indice switch.
# LAMMPS convention is NOT the Voigt notation.
variable aves1 equal "ave(f_VEC[1])"
variable aves2 equal "ave(f_VEC[2])"
variable aves3 equal "ave(f_VEC[3])"
variable aves4 equal "ave(f_VEC[6])"
variable aves5 equal "ave(f_VEC[5])"
variable aves6 equal "ave(f_VEC[4])"

# Computing the covariance through the <s_{i}s_{j}>-<s_i><s_j>
# is numerically instable. Here we go through the <(s-<s>)^2>
# definition.

# Computing difference relative to average values
variable ds1 vector "f_VEC[1]-v_aves1"
variable ds2 vector "f_VEC[2]-v_aves2"
variable ds3 vector "f_VEC[3]-v_aves3"
variable ds4 vector "f_VEC[4]-v_aves4"
variable ds5 vector "f_VEC[5]-v_aves5"
variable ds6 vector "f_VEC[6]-v_aves6"

# Squaring and averaging
variable dds1 vector "v_ds1*v_ds1"
variable dds2 vector "v_ds2*v_ds2"
variable dds3 vector "v_ds3*v_ds3"
variable vars1 equal "ave(v_dds1)"
variable vars2 equal "ave(v_dds2)"
variable vars3 equal "ave(v_dds3)"
variable C11 equal "v_Pa2GPa*(v_C1*f_CB[1] - v_C2*v_vars1 + 2*v_Rhokbt)"
variable C22 equal "v_Pa2GPa*(v_C1*f_CB[2] - v_C2*v_vars2 + 2*v_Rhokbt)"
variable C33 equal "v_Pa2GPa*(v_C1*f_CB[3] - v_C2*v_vars3 + 2*v_Rhokbt)"

variable dds12 vector "v_ds1*v_ds2"
variable dds13 vector "v_ds1*v_ds3"
variable dds23 vector "v_ds2*v_ds3"
variable vars12 equal "ave(v_dds12)"
variable vars13 equal "ave(v_dds13)"
variable vars23 equal "ave(v_dds23)"
variable C12 equal "v_Pa2GPa*(v_C1*f_CB[7] - v_C2*v_vars12)"
variable C13 equal "v_Pa2GPa*(v_C1*f_CB[8] - v_C2*v_vars13)"
variable C23 equal "v_Pa2GPa*(v_C1*f_CB[12] - v_C2*v_vars23)"

variable dds4 vector "v_ds4*v_ds4"
variable dds5 vector "v_ds5*v_ds5"
variable dds6 vector "v_ds6*v_ds6"
variable vars4 equal "ave(v_dds4)"
variable vars5 equal "ave(v_dds5)"
variable vars6 equal "ave(v_dds6)"
variable C44 equal "v_Pa2GPa*(v_C1*f_CB[4] - v_C2*v_vars4 + v_Rhokbt)"
variable C55 equal "v_Pa2GPa*(v_C1*f_CB[5] - v_C2*v_vars5 + v_Rhokbt)"
variable C66 equal "v_Pa2GPa*(v_C1*f_CB[6] - v_C2*v_vars6 + v_Rhokbt)"

variable aC11 equal "(v_C11 + v_C22 + v_C33)/3."
variable aC12 equal "(v_C12 + v_C13 + v_C23)/3."
variable aC44 equal "(v_C44 + v_C55 + v_C66)/3."

print """
C11 = ${aC11}
C12 = ${aC12}
C44 = ${aC44}
"""

C11 = 2.43259565471636
C12 = 1.31707509524277
C44 = 1.46049163278052

Total wall time: 0:00:28
